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SUMMARY 


This report is a compilation of the research in orbit determination and 
optimization in space trajectories carried out by AMA, Inc, , under contract to 
Lyndon B. Johnson Space Center, covering the period February 1972 - September 
1975, 
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INTRODUCTION 


Analytical Mechanics Associates, Inc., under contract to the Lyndon B. 

Johnson Space Center, acted in the capacity of consultants in the areas of orbit 
determination, optimization techniques and trajectory design for manned space 
flights. In this capacity, several reports were generated and are included in the 
text of this final report. 

A brief description of each report is included here. 

(1) Multi-Spectral Scanners 

This report contains a recommended method applying optimization techniques 
for estimating the most likely source stimulus for a given sequence of multi-spectral 
signals obtained from a flight scanner passing over an area. The report contains well 
known statistical discrimination techniques for estimating and refining information 
from multi-spectral scanners. 

(2) Non Singular Earth Gravity Acceleration for Space Shuttle 

This report contains a computer routine for computing the earth gravity accelera- 
tion designed for space borne computers. The method is both time and computer core 
efficient. 

(3) Drag Acceleration as a Navigation Aid 

This report develops a technique for obtaining a more accurate estimate of the 
vehicle state during reentry blackout from a drag acceleration measurement by in- 
cluding an adaptive atmospheric temperature lapse rate as an additional model para- 
meter. 

(4) Roll-Modulated Lifting Entry Optimization 

This report develops the optimal technique for a roll-modulated lifting reentry. 
Thus, it provides valuable insight for obtaining practical reentry guidance laws using 
roll modulation. 
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(5) Minimum- Variance Linear Estimator for Non-Linear Measurements 

This report derives the minimum variance estimator for modifying the Kalman 
filter using range-rate measurements during highly non-linear flight regimes, 

(6) Approximate Matrix Inverses 

This report derives a rapid method for obtaining an approximate matrix in- 
verse for use in flight computers when a full inverse is too time consuming and not 
numerically compelling, 

(7) A Variable-Metric Algorithm Employing Linear and Quadratic Penaltie s 

This report develops a variable metric optimization algorithm for accelerated 
search for optimization problems with linear and non-linear constraints. 
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AN APPROACH TO VIEWING OP 


MULTI-SPECTRAL SCANNER DATA 

Consider the problem of presenting mnlti-spectral data to a viewer via a 
cathode-ray tube, say a conventional TV or, perhaps, a color TV. Spectral dis- 
crimination is to be used in some way to aid distinction between objects of main 
interest (targets) and other objects (background). Exploitation of the human T s 
pattern-recognition capability is, of course, the main attraction of the viewing ap- 
proach, Preprocessing of the data as well as the data for ground-truth tracts is 
assumed on the basis of purely Spectral discrimination computations, from which 
at least approximate models of target and background spectra and probability 
density distributions are available. 

With pt , , 11 the intensities of signal in the n measurement fre- 

quency bands and the probability densities of target and background denoted 

f T^ ,U l ’ 3 ^ ’ » , as P er 1 » *&© rati0 of tlie two 



is the criterion often employed for purely spectral discrimination. Thus, if , 
the sample is interpreted as a target signal. The threshold C is set to admit a 
specified percentage of known targets. The approach contemplated is preprocessing, 
perhaps approximate, to obtain models for and including numerical values 
of means and covariances of components, then evaluation of the likelihood ratio K, 
for each data sample and use of it as a signal for a display option. 
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Mean, and covariance for the target probability density distribution are com- 
puted from ground truth data, and in many cases a Gaussian representation 'Will 
suffice. The background sometimes may be approximated by a Gaussian distribu- 
tion, but more commonly 'will better by taken as the sum of two or three or several 
such, each term representing some prominent ingredient. To get means and co- 
variances for the background, the total signal must be processed and those parts 
clearly target contribution screened out by a coarse criterion such as mean = 2c. 
With this deletion, the mean and covariance of the remainder can be computed and 
the distribution represented as a sum of Gaussian contributions. 

Feeding the likelihood ratio to a video display would light up the target areas 
and leave others dark. If there are two or three targets, a color display could be 
used. For the purpose of estimating the likelihood ratio for a particular target, the 
other targets are treated as components of the background. The display should, of 
course, also function normally, with choice of displaying the picture in any of the 
bandwidths on option, or perhaps there should be a second display for simultaneous 
viewing of raw data. In this fashion, the likelihood ratio would not be used for classi- 
fication according to a black-or-white threshold criterion, but would furnish shades 
of grey (or green) for display and the viewer, thus assisted, would do the classifying. 
It is difficult to anticipate whether the approach sketched here might find its best ap- 
plication in preliminary editing of massive amounts of data, perhaps using off-the- 
shelf density models from previous reductions, or whether it might instead excel for 
intensive efforts on infrequent difficult cases. 
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INTRODUCTION 


Tlie Shuttle computer requires a nonsingular gravity acceleration routine 
for polar orbits. The folio-wing formulation, based on Reference 1, is carried 
out in detail for an Earth model consisting of , 3 ^ , , and & . 
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EQUATIONS 


■where 


The nonsingular potential in Earfch-fixecl coordinates is given by 

* ' “ V " 1(f) 2 J 2 A 2, 0« " |(fj J 3 A 3, GW - ? (f f J 4 A 4,0< U > 
4 £-(*) Vt0 2!¥ ,t) 

(1 = central body gravity constant 

a = Earth radius 

, 2 2 2. 4 
r = (x + y + s ) a 

x 

S = — 
r 

t = y - 
r 

z • 

u = — 

r 


, 1 d 1 * 3 ,2'i 

a. . = — — — -rrr (ii - 1) 

1,3 2 l i! du l 3 
R 2 (s,t) = s 2 ~t 2 
I 0 (s,t) > 2 s t 


( 1 ) 
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The acceleration vector in the body- fixed system is given by 
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-( 
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where 1, j , k are unit vectors in the x, y, a directions, respectively. 


For the specific Earth oblateness coefficients J , , J ( . G , and 

2 t> 22 

^22 5 ^ 1C P ar ^^ s are gi-^en below. 
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A simple computer code is given in the next section. 


N ON SINGULAR EARTH OBLATENESS ROUTINE 

OBLAT(R,ACC) 

DIMENSION R(3), ACC (3) 

COMMON EMU, a, 32, 33, J4, C22, S22 
COMMENT EMU - central mass coefficient 
COMMENT a - Earth radius 

R2 > R(1)*R(1) 4- R(2)*R(2) + R(3)*R(3) 

R = SQRT(R2) 

EINV = 1.0/R 

B = a *RINV 
C = EMU/R2 * B 
B2 = B *C 
B3 = B*B2 
B4 = B*B3 
S = R(1)*RINV 
T = R{2) *RINV 

U = R(3) *RINV 

AS = B2*6.0*(S*C22 + T*S22} 

AT = B2 * 6. 0 * (S * S22 - T * C22) 

XJ2 = U*U 

U3 = U*U2 

U4 U*U3 

A2 = (3.0*U2 - 1.0)/2.0 

A21 = 3. 0 He XT 

A3 - (5. 0 *U3 - A2l)/2. 0 

A4 = (35. 0 * U4 - 80. 0 * U2 + 3, 0)/8. 0 


A31 = 

A41 = 

AU 
AR 

AR 

ACC(l) 

ACC(2) 

ACC(3) 

END 


(15,0*U2 - 3. 0)/2.0 
(35. 0 *U3 - 15,0*U)/2.0 
- B2 * J2 * A21 - B3 * J3 * A3 1 - B4 * J4 * A41 
3.0*B2*J2*A2 + 4.0+B3*J3*A3 + 5. 0 *B4 * J4 * A4 
- B2 * 3. 0 * (C22 * (S *S - T * T) -I- S22 * 2. 0 * S * T} 

AR - S* AS - T * AT - U*AU 
AR*S +AS*RINV 
AR * T + AT * RINV 
AR * U AU^RINV 
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NOTATION 


p atmospheric pressure 

p atmospheric mass density 

T atmospheric temperature (degrees Kelvin) 

m mean molecular weight 

\v 

q universal gas constant 

h altitude 

B vehicle position vector 

• 

B vehicle velocity vector 

V_, vehicle velocity vector relative to the atmosphere 

g acceleration of gravity (assumed constant in the atmosphere) 

AV delta velocity readout vector of the accelerometers 

At accelerometer readout count time interval 

A effective vehicle area for drag 

m vehicle mass 

CO Earth angular rotation rate 

M Mach number 

c speed of sound in atmosphere 

r_ mean Earth radius 

E 

e eccentricity of Earth 


ii 


0 ' 


geocentric latitude 


0 geodetic latitude 

z polar component of vehicle position vector 

y ratio of specific heats 

r magnitude of radius vector 

v R magnitude of relative vehicle velocity 


iii 


INTRODUCTION 


The drag acceleration is presently being considered as an observation type 
for navigation as a substitute for an altimeter reading during radio blackout in re- 
entry. The measurement is somewhat degraded in the event the atmospheric 
density or the aerodynamic drag characteristics become uncertain. This report 
resolves the density variability by using a model of the atmospheric density from 
SO Ion to 32 km in the form of a layered atmosphere characterized by a sequence 
of altitude intervals with piece-wise constant temperature lapse rates. The intent 
is to use this parameterization to enable the navigation filter to determine the 
vehicle state as well as the atmospheric density variation. In this manner, a 
more accurate determination of the state and covariance cross-correlation will 
exist at blackout termination when more effective observations can again be made. 
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ATMOSPHERIC MODELS 


Below 100 km , the Earth's atmosphere seems to obey the perfect gas 
law for a constant mean molecular weight. 


m p 
w 

P ~ qT 


where 


m 

= constant 

q 


(i) 


Under the assumption of a temperature variation which is piece-wise linear with 
altitude, we have 


h £ h ^ h 
1 2 


T = T^) + T^h-h^ 


( 2 ) 


Since the equilibrium pressure variation is given by 


dp 

S = 


(S) 


the density between h and h is given by 


T*(h ) 

P = P( h i) [} * T(h x ) 


-C[gm w /qT l (h 1 )] + l} 


(4) 


Since the equilibrium atmosphere experiences temperature gradient re- 
versals, there exist altitude regions of 5 km or so during which constant 
temperature persists. In such intervals, we have 


Hr- 2 


T(^) = T = T(h 2 ) 

-Cm - v g/TT(h DlCh-hJ, 
p = p(h x ) e 

The 1962 U. S. Standard Atmosphere (Eef. 1) lists the following nominal 
values of the temperature lapse rates 


TABLE I 


Altitude 

Ion 

Temperature 

^Kelvin 

Lapse Bate 
°Kelvin/km 

Density 

kg/m 3 

32 

228.65 

2.8 

1,3225x10" 

47 

270. 65 

0 

1.4275x10" 

52 

270.65 

-2.0 

7.5943x10" 

61 

252, 65 

-4.0 

2. 5109x10" 

79 

180.65 


2.001 xlO" 


( 5 ) 


Using the above four-layer model, it is possible to construct an error model for 
density variations and require that the filter solve for the parameters during the 
descent blackout. 
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DRAG ACCELERATION 

The observation of drag acceleration over the short time interval during 
•which the accelerometers read out the delta v count (0. 5 seconds) is given by 


, i 2 _ A 
D.A. = 2 P V R C Dm 


av-v 


R 


% At 


(S) 


While the representation of the observation appears to provide the ability to 

estimate both position and velocity, the uncertainty in modelling the aerodynamic 

coefficient C as a function of Mach number and angle of attack make this ap- 

2 A 

proaeh undesirable. We will compute the quantity % v^ ~ and produce the 
pseudo-observation of density 


AV'W 


p = 


R 2m 


v| At C D A 

}x 


(7) 


From the knowledge of the vehicle state, we compute 


V = R-^xR 
R 

s 


R 

Cl 

°D 

M 


<V V E> 


= k 10 


(Earth 1 s angular rotation vector) 


(^(M , a) (drag coefficient as function of 
Mach number of angle of attack) 

\ 

c 

[r iw]* 


( 8 ) 


m 

w 


T(h) = 


a = 


T(h x ) + T’/h^th-h^ 


cos 


„ • V 
-1 x R 


V. 


(angle of attack) 


R 
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Following the previous section, we model the density as a function of the 
geodetic altitude, h. The altitude is the height above the vehicle sub-satellite 
point. The Earth is modelled as an oblate spheroid with eccentricity, e. From 
Ref. 2, the altitude is given by 


i 2 2' 

h = r cos(0-0') - r (1-e sin 0) 


( 9 ) 


The geocentric latitude is defined as 


* jl „ Z 

Sill 0 = — 

r 


The geodetic latitude, accurate to twelve digits, is given by 

0=0'+ a^(r 7 e ) sin 2 0 r + a^(r , e) sin 40’ + a g {r , e) sin 6 0 1 


+ a g (r , e) sin 8 0* 


where 


\ = T l0M (512e 2 + 12Se 4 -HS0e 6 + S 5e 8 ) + (^)|j( e 8 + e S ) 


VTS (54e 4 4«e 6 + W. 8 ) + (3) 4 e 4 + 2 B 6 + e 3 ) + ^(2) 


2 3 

.8. . f\\ 1 „ 4. „ 6 . 8, . 15 e / r E' 


8 


-(?) 


£ 

16~ 


2 r 3 

r E 3 ’ . 6 , _ 8, 3 / ,3 E\ . 6 , 8 35 /TEN ..8 8. 

a e = T Io24 (4e +5e (e * e ) + t£S\ 77 (4e ' I ' 3e } 


(9a) 


(9b) 


(9c) 


a 8 2048 


[-«T + m (t) -“*(?)**”(?)] 




,2 E' 


in- 5 


2 2 
e - 2 £ - £ 


£ 


1 

297.3 


The above equations are sufficient to enable one to estimate the density. 
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PARTIAL DERIVATIVES OF THE DENSITY OBSERVATION 


We take as our state 


R 

+ 

R 


P(*\) 

T'(K) 

T(K) 


vehicle position vector 
vehicle velocity vector 


th 


density at lower altitude of the i layer 


th 

temperature lapse rate of the i layer 


.th 


temperature at the lower altitude of the i layer 


The partials of the density with respect to each component of the state 
are given by 

x * * bp( h.) . 5T’(h.) 

bp _ bg_ Mi , bp r bp i 

bp bh bjB ‘ r bp(h.) bp bT r (h.) bp 
bp «<*i> 


bT(li.) bfi 


( 10 ) 


The partials of the density with respect to the altitude, h, are given by 
For T f (h ) f 0 




T'ft) 


.q T (h.) 


T(h.) 


T(hJ 

(Ha) 


For T'(Ii} = 0 


>, S m 

&jo = & w 

bh p q T(h.) 


(Ub) 
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The partial of the density with respect to p (h.) is given by 

bp = P 

bpiK) p(h.) 

The partial of the density with respect to T f (h.) is given by 
For T 1 (h.) ^ 0 

h~h. 

W 

For T^h.) = 0 
bT r (h.) 


bp 

bT'(h.) 


pjjtn 


{(i + 


T'(h.KIi-h.) 


T(h.) 


ty-gr:)- 


1 + 


gm 

w 

<1 T r (fl.) 


3 T ,? '(a.) 


1 + 


T'Ch.) 

W 


The partial of the density with respect to T(lr) is given by 


For T 1 (h.) ^ 0 


gm 


w 


OT(h.) T'{h.) 


T (h*)(h-h.) _ w 


- CIg bi /q T t (h.) ] + 2} ^ 


-I— 


T^) 


For T* (h.) = 0 


* g m 

bp _ w 


bTfh.) 2,, . 

i q T (h.) 
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(13a) 


(13b) 


(h.)(h-h.) 

T 2 {h.) 

(14a) 


(14b) 


The partial of the altitude, h , with respect to the state is computed on 
the assumption that 0=0*, and is given by 


bh _ 

r 

Ti 

2.2' - 
e sin 0 

12 



bR 

L 1 " r 

(1-e 2 sin 2 0) a “ 

1 r + 

r 

(1 

bh 

bh 

bh 

bh 


n 

O' 1 

i 

bT’ (h.) 

bT(ll.) 

^P(^) 


U 


2 . 2 , 
e sin 0 

2 


i h 

.a 


( 15 ) 


The partial of p(h.) with respect to the state is 

bp(h.) bp (h.) bp(h.) bp(h.) 

= = ^(b.) = bT(h.)' = ° 

bpch.) 

Sp(h.) 1 


(16) 


The partial of the temperature lapse rate with respect to the state is 
given by 


bT T (h.) bT’(h.) bT'(K) bT l (h.) 

bK " _ “ bp (h.j" ~ bT(h.) 

bT'Oi.) 

bT’(h.) 1 


(17) 
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4 


The partial of the temperature at h. , T(h.) , with respect, to the state is 
given by 


bT(h.) 

bit" 


bT(h.) 

bR 


bT(h.) bT(h.) 

bp(Xj = bT l (h.) 


bT(h.) 

bT(h.) 


(IS) 


This completes the partials. 

■Whenever [T T (hj)J^ 0.5 /km, set T (h.) — 0. 
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NAVIGATION FILTER UPDATE EQUATIONS 


Reference 3 contains the update equations for an onboard navigation filter 
for an 18-element vehicle state vector consisting of the vehicle position vector, 

R , the vehicle velocity vector, R , the vehicle gyro tilt error vector, 6 , the 
gyro tilt rate error vector, 0 , the accelerometer scale factor error- vector, k , 
and the accelerometer bias error vector, b . To this state we add the scalars 
p (h.) , T 1 (1\) , and T(h.) . Thus the new state is a 21~element state error vec- 
tor, X. Following Ref, 3, the update equations following each pseudo-density 
observation are given by 

T "1 

X + = X + C P T (PCP T +Q) (p.-p ) (19) 

v ' T>b 'comp' v 


where C is the 21x21 covariance matrix of the errors in the estimate of the 
state vector, P is the 21x1 vector of the partials of the pseudo-observation 
with respect to the state X, and Q is the observation noise. We have 


•p - F M &P bp b_p bp bp 

* L bR 5 1 be ’ 1 bk ’ bb ’ bp(h ) * &T'(h ) 


1 

bT(h.) J 


( 20 ) 


m 

The partials of the density with respect to E, R, p (lr) , T T (h.) , and 
T(h.) are given in the previous section, and the remaining partials are all zero. 

The residual is given by 

AV • V 2 m 

Ap - - pCh, p(h), T(h), T'(h)] (21) 

O _ . L L L 

\ c d a 

The value of Q recommended for simulation study is 

Q = *01p 2 (22) 
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RoU-Modulated Lifting Entry Optimization 


ORIGINAL PAGE IS 
QE POOR QUALM 


Henry X Kelley* 

Analytical Mechanics Associates Inc., Jericho, N.Y. 

AND 

Henry C. Sullivan! 

Lyndon B. Johnson Space Center, Houston, Texas 


The equations of lifting entry are examined for fixed angle-of-attack vehicular motion with path control via 
rail modulation of lift. A complication arising with this is nonconvexlty of the heliograph figure, which mnkes the 
application of standard variational techniques inadvisable unless the problem is first relaxed, i.e,, a related problem 
is defined with a hodograpli figure that is the convex hull of the original. This leads to a hew system in new 
variables that is apparently innocuous in its simplicity; the linear elements of the convex hull, however, are associated 
with singular extremal subarcs and their attendant difficulties. The singular extremal for minimum-heating 
symmetric flight with final time and downrange open is simple. Two order-reduction approximations are considered, 
which may include intervals of two-dimensional motion as subarcs. One.of these approximations relegates turning tD 
initial and terminal boundary-layer maneuvers; the other is analogous to the aircraft energy-maneuvering model. 
Some computations for a space shuttle orblter configuration are presented. 


Nomenclature 

D - drag 

E = specific energy 
g 0 - acceleration of gravity 
H = variational Hamiltonian 
L = lift 

0 = total heat load 
6 = heat rate 

r - radius 

r p = radius or the Earth 
V — velocity 
IV= weight 

y = flight path angle to horizontal 
A = longitude 
7. = Lagrange multiplier 
/r = bank angle 

; = relaxation interpolation variable 
a = relaxation control variable 
g = latitude 

1 = heading angle to south 


State Equations 

W ITH r radius, y path angle to horizontal, E — (V 2 /2i] 0 )- 
(r a 2 /r) specific energy, / heading angle to south, <f> latitude, 


A longitude, and p bank angle, the equations of state are 

r = V sin y (1) 

£ = —DV/W (2) 

y = (g 0 Lco5 pjWVj—lgy 2 jVr 1 ) cos y + (7/r) cosy (3) 

'/. ~ {g 0 Lsu\ [tjWV cosy)— (F/r) cosy sin y.tan 0 (4) 

r= - (H/r} cos x cos y (5) 

A = ( V sin /_ cos y/rcos <j>) (6) 

£=Q{E.r) ( 7 ) 


Presented as Paper 72-933 at the AlAA/AAS Astrodynamics 
Specialist Conference, Palo Alto. Calif., September 11-12, t972; 
submitted October 6. 1972; revision received March 13, 1973. 
Research supported in part by Lyndon 13. Johnson Space Center under 
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Index categories: Entry Vehicle Mission Studies and Plight 
Mechanics; Navigation. Control, and Guidance Theory, 

* Vice President. Associate Fellow AIAA. 
t Aerospace Technologist. Member AIAA, 


The first six equations are particle-dynamics equations of motion 
for coordinated maneuvering (zero side-force). The last equation 
is the total heating integral Q in differential Form. Lift L and drag 
D are functions of E and r only; angle of attack is assumed 
constant. {IF trim were to vary with the Mach number, the aiHe 
of attack would itself be a function of E and r.) Inequality con- 
straints on dynamic pressure, normal load factor, and local 
temperatures are in the problem statement. 

Roll Modulation 

Entry at essentially constant angleof attack has been employed 
for such vehicles as the Apollo Command Module, with 
consequent simplification of longitudinal control. The desired 
vertical component of lift and a desired average out-of-plane 
component are obtained by bank reversals, square-wave fashion. 
In the particle-dynamics model, this includes the theoretical 
possibility of “chattering,” since rigid-bo dy rolling dynamics have 
been neglected. There is design interest in roll modulation for 
advanced vehicles such as the Earth-orbital shuttle, even though 
a longitudinal control system will be featured, since design 
compromises may force a narrow range of trim angle of attack. 
Tlius, constant angle-of-attack operation is of interest as a 
limiting case for the shuttle entry problem. 

Control Relaxation 

In the version of the problem with angle of attack controllable 
within bounds; the figure in hodograph space (£, y, j?) that is 
traced out by varying the controls a and p over their complete 
range{ContensouV l Domain ofManeuverability”) 1 isnotconvex. 
Operation at points within the figure, which is a paraboloid for lift 
linear and drag quadratic in i, can be approximated by chattering 
control operation, square-wave fashion, but cannot actually be 
attained with piecewise continuous controls. In such circum- 
stances, it is usual to consider instead a related problem with 
different control variables that attain the convex hull of the 
h odograph figure ; this is thc'Vclaxcd" problem. 1 2 The relaxation 
for the variable angle-of-attack case is sketched in Ref. 3, In the 
present case of fixed angle of attack, the figure is an ellipse. 
Relaxation makes the disk within this ellipse attainable. 

Relaxation may be accomplished for a general state system of 
the form 

& -/(.v.ii.t) (8) 
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Fig. 1 Minimum-heating trajectory in altitude/velocity chart. 


by replacing the system by 

x —fix, n, . i) + £[/(*■ tt z , t) —f (x. tt,, f)] (9) 

in which the right members are linearly interpolated between 
values for control ir t and control u z . Here £. 0 £ £ g l, is an 
interpolation parameter. The control variables of the relaxed 
system are the vectors it, and u 2 and the scalar £. In the present 
application, the desired goal of attaining the interior of the ellipse 
can be accomplished with fewer variables, namely by introducing 
an additional control variable c, 0 g a g 1, multiplicative on 
L in the y and y state equations 

V = (&,L<7 cos lt/WV)—ic] 9 r a 2 /Vr z ) cos y-r(Vfr) cosy (3a) 
y = [g a U t sin jijV/V cosy)- (Vjr ) cosy sin y tan r}> (4a) 

Singular Arcs of the Relaxed Problem 

The appearance of the control variable c linearly in the right 
members of the state equations indicates the possibility of si ngular 
res in the solution of optimal entry control problems. This 
possibility may be investigated by formation of the usual 
Hamiltonian //, setting cH/cg — 0, and pursuing the 
consequences. 

dH/ca — (cj 0 L/W V)[/. ? cas /r+ a, (sin /i/cosy)] = 0 (10) 

dH/dfi - (g e Lt?/WV)[— /. y sin (cos p/cosy)] = 0 (11) 

Left members of Eqs. (10) and (11) must vanish independently. 
Since these are linearly independent, it follows that both and 
}. y are zero along the arc. 

The system is already in the canonical Form of Ref. 4; thus the 
variables y and y are control-tike along singular arcs. A similar 
result could have been obtained by noting that a sin n and 
g cos could be taken as new control variables in the neighbor- 
hood of a singular arc ForO < a < 1. Desired variations in y and 
■/ can be realized by varying these, as long as the magnitude of the 
desired variations is sufficiently small as not to encounter 
saturation of the a bounds. 

With y and y regarded as controls, the problem simplifies to 
flight in the plane of a great circle. Without loss of generality, take 
g — ijj = O.andy = jt/ 2 for study of this two-dimensional motion. 


and the state equations became 

r=Fsiny (12) 

&=-DV; W (13) 

A = V cosy/p (14) 

0 = ()(£.r) (15) 


In the special case of downrange open (final A unspecified for 
initially equatorial flight), the control variable y enters only 
Eq (12) and the variable r becomes control-like along singular 
arcs as the form with Eq. (12) deleted is again canonical. If final 


time is open, there is analytical advantage in casting E in the role 
of independent variable; furthermore, the steady decrease of E 
makes this interchange feasible for entry applications. 

dQ/dE =. - WQ/DV (16) 

The singular extremal is defined by stationary points of the right 
member of Eq. ( 1 6) regarded as a function of r at various E levels. 
The generalized Legendre-Clebsch 4 " test requires that the 
stationary value of the right member of Eq. (1 6) as a function of r 
be a maximum and Q/DV minimum. 

Reduced Relaxed Problems 

The relaxed problem presents computational difficulties 
because of singular arcs; approximations are therefore of mrv-e 
than usual interest. Possibilities offered by singular perturba- 
tion procedures 5-7 are discussed in the following paragraphs. 
If nearly symmetric flight were assumed, a singular perturba- 
tion approach designating latitude, longitude, and heating 
as variables of a reduced solution (i.e„ solution of a reduced- 
order approximation problem) would seem attractive. This would 
relegate turning and altitude transitions to cc rrective boundary- 
layer transients near initial and terminal points. Energy is 
chosen as the independent variable. The reduced problem is of 
the great-circle type. 

The great-circle reduced-order system for the approximation 
that combines altitude and heading transients takes the form 

dA/dE = - W/Dr (17) 

dQ/dE = -WQ/DV (18) 

The order-reduction procedure used is the same one examined 
and employed in Ref. 6. An upper bound on the control variable 
r of the reduced problem is furnished by the control bound a 
of the original problem together with Eq. (3a) and y = y = 0; 
a lower bound is provided by state inequalities on panel 
temperatures and acceleration loads, handled in penalty function 
approximation in the computations next described. Use of the 
model given by Eqs. (17) and (18) is limited to problems for which 
downrange is specified as greater than, or equal to, the down- 
range-open value; for smaller specified downrange, the singular 
extremal fails the generalized Legeudre-Clebsch test and a zigzag 
competitor is optimal. 

A less drastic approximation using singular perturbations 
would treat heading as well as latitude, longitude, and heating 
in a reduced problem. This would idealize only the altitude 
transients as fast (with respect to energy change) compared to the 
other transitions. It is the same as aircraft energy approxima- 
tion. 3 - 6 Energy approximations have previously been examined 
for atmospheric entry of a variable anglc-of-attack vehicle. 7 
No complications arising from the relaxed model are anticipated 
using this approach. Evidently a solution for the reduced-order 
fixed angle-of-atlack problem consists generally of a turning arc, 
a great-circle time-open arc, and, if final heading is specified, a 
terminal turning arc. 

Computational Results 

Data for a delta- wing space-shuttle orbiter configuration were 
used for some sample computations with the model of Eqs. (17) 
and (18). The angle of attack was fixed at 30’. Inequality con- 
straints on normal load factor and numerous panel temperatures 
were incorporated by using penalty functions. 

A minimum of the Hamiltonian consisting of a linear com- 
bination of the right members of Eqs, (17) and (18) plus penalties 
was found by one-dimensional search. With downrange open, the 
minimum always occurred at the lower bound on altitude 
furnished by the toad factor and temperature constraints (see Fig. 
1 ). With downrange specified at values exceeding the open value, 
the minimizing altitude was found to be the upper bound value 
(?r = l)during the latter part ofthe trajectory (see Fig. 2). As range 
requirements were increased, numerical results indicated the 
possibility of more than one switch between altitude bounds. 
The Hamiltonian function for the down range-specified case of 
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Fig, 3 Hamiltonian vs altitude at several energv levels far dawnmnge- 
veiocitt magnitude • ft*s specified case. 


Fig. 2 Minimum-heating downrange-specified trajectories in altitude/ 
velocity chart. 


Fig. 2 is sketched vs altitude in Fig. 3 for several energy values. 
The sign of the second derivative H rr would seem to indicate 
nonconvexity and a need for further relaxation. However, one 
recalls that rin the roleof control variable is not the real thing but 
the result of an order-reduction approximation amounting to 
assumed instantaneous vertical dynamics. This implies that weak 
as well as strong minima should be considered; hence, that 
transitions determined according to absolute minimum ff, as in 
Fig. 2, arc somewhat arbitrary. 0 Boundary-layer transition 
fairings at discontinuities in r, as in Ref. 6, arc needed for 
consistency in degree of approximation of die control, but they 
contribute nothing to the performance index in this 
approximation. 

When heating was heavily weighted compared to down- 
ranging, values of <7 were found to be below unity indicating a 
need Tor roll modulation in two-dimensional flight. However, 
solutions with downrange heavily weighted ride the upper bound 
o’ = t at low energies and at near-orbital energies, 

Results obtained by a conjugate gradient method that used a 
particle-dynamics model are shown for comparison. The cross- 
range was specified at a somewhat challenging value of about 
1300 nm. The conjugate gradient formulation did not employ a 
relaxed model and was unsuitable for nearly symmetric flight 
cases. It exhibited poor convergence that was, perhaps, 
attributable to the absence of convexity. Nonetheless, the result of 
Fig. 2 seems of interest for the qualitative similarity of the 
altitude history with the great-circle model. This was obtained 
using as a first guess a trajectory which had been forced to follow 
the lower bound representing temperature limit approximately. 
The comparisons suggest that the idealization of early heading 
and altitude transitions followed by altitude control based mainly 


on heating and down-ranging may warrant further investigation. 
A separate treatment of the initial transition as a boundary layer 
in which altitude, path angle, and heading motions are fiist (with 
respect to E changes) could be carried out along the lines of that 
for aircraft altitude transitions in Ref. 6. 

Conclusion 

Attention has been directed to relaxation and its consequences 
for the fixed unglc-of-altack atmospheric entry problem. Two 
reduced-order approximations for the resulting system of 
equations have been briefly examined and appear to warrant 
additional study. 
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Let be the best estimate of the vector state, and x the true state vec- 
tor. Then the vector error in the best estimate ls given by 

e = x - x (1) 

oo 

Let y be a scalar measurement which is a nonlinear (quadratic) function 
of the vector state x, contaminated by white noise. Then, the true measure- 
ment is given by 

a T T 

y( 30, rj) = y(x } +y e + % e y e + 77 (2) 

' " O X O O XX 0 

We seek a minimum variance linear estimate of x of the form 

x-x q = K[y(x,??) -y(x)] (3) 

where K is the linear vector gain. The expected change in the error in the 
estimate is given by 


e = e Q ~ 


(4) 


The variance of e is given by 


E(e,e T ) = E(e Q ,e o T ) -KE(Ay,e o T ) -E(e o ,Ay T )K T 


+ K E(Ay,Ay T )K T 


(5) 


The minimum variance estimate of e over all linear K gain vectors is given by 


K 


^(Ay , ej ) 


mv . T. 

E(Ay , Ay ) 


( 6 ) 


V-l 


X 


If we now make the assumption that 


and that 


T T 

E(e , 5 ]) = E( e , e y e ) = E(-n ,e y e ) = 0 
' o " v o o J xx o' K 1 o J xx o' 


E{e ,e T ) = P 
o o 


EC??,?? ) - Q 


we have 




K = 


x 


T ITT 

y Py + Q+TE[e y e,e y e) 
x x 4 v o xx o o xx o' 


Prom matrix algebra, we have, for any vector <£, 


£ T y [It*’) y l trace [y (££ T )y (£-i- T )] 

xx xx xx xx 


We seek the expected value of the trace of the square of the matrix A , wher 


A = y (e e ) 
xx' o o ' 


Given any nxn matrix A, the trace of its square is given by 

N N 

trace A 2 = V t a., a.. 

ft ft l 3 J l 

i=l j-1 

Every element of the matrix A is given by 


N 


a ij = 1 hi e 0k e 0j 


k=l 


til 

where y is the ilc element of y . 
ik xx 


We seek the expected value of 


N N N N 

E Z l I I y tk y U B 0k e 0i ®0-t, e 01 


1=1 j=l b=l £=1 


( 13 ) 


For a normally-distributed random variable with zero mean, we have 


Ete 0k e 0j e 0l e 0i* ECe Ok e Oj^ E(e 0^ e Oi^ +E<e Ok e O^ :E(e Oj e 0i* 


+ E(e 0k e Oi )E{8 03 8 <U> 


(14) 


Let 


p„ = E ( e e ) 
i] Ol O 3 


(15) 


Then Eq. (13) becomes 

N N N N 


III I y ik y j 

i=l j=l k=l 1=1 . 


-L ^ ^kj P li ‘ P k-t *ji + P ki ^3 


(16) 


If we define the matrix C to be the expected value of the matrix A 


C = E(y e e ) = y P 
xx 00 xx 


(17) 


Equation (17) may be written as 

N N 

y y c.. c.. +c.. c.. +c.. c..) 

A A U ,U u n u 33 


1=1 j=l 


( 18 ) 


from which we obtain 


AQ = ~ E trace (A 2 ) 


1 2 1 2 
— trace (C )+— (trace C) 

Ik 


V-3 


(19) 


For the purposes of computation, we recommend that we first compute 

the nonzero elements of C=y P and then form AQ. 

xx 

N N N N 2 * 

AQ = f I c u + 1 I SjVitfiv] < so > 

i-1 i-1 j =1-5-1 l-l 


Example 


Let y depend only upon position (e. g. , in a range measurement, in 


DML, or in angle measurement) . Then y is an upper 3x3 and the trace 

12 xx 

of the expected value of — A is given by 


Let 


= c (3^3) = 


,1 

'4 


1 . „2 


11 

C 12 

C 13 

21 

°22 

C 23 

31 

C 22 

C 33 . 

+C 22 

+ C 33> 

+ < C 12 < 


( 21 ) 


12 21 13 31 23 32 - 


+ 4 (C 11 + C 22 + C 33 J 


( 22 ) 
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SUMMARY 


This report derives a procedure for a rapid determination of an approxi- 
mate matrix inverse for use on real-time on-hoard computer control systems. 
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INTRODUCTION 


On-boarcl computers often require matrix inversions during real-time 
computer control system computations when the cycle time does not allow for 
precise solution. The procedure outlined here yields an approximate solution 
which can be executed in considerably less computer time. 


DERIVATION OF THE APPROXIMATION 


Let A be an nxn matrix, and x and 
an approximate solution of the problem: Given 


y be nxl 
A and y , 


vectors. We seek 
find x , when 


Ax " v 


( 1 ) 


Let B be the matrix formed of the diagonal elements of A. 


b„ = 6,. a., 
ij -3 13 


(2a) 


Then 


A = B B 1 A 


(2b) 


The solution of Eq. (1) is given by 


~1 1 -1 
x = (B A) B y 


( 3 ) 


Let 

C = b” 1 A 

and (4) 

z‘ - B y 

The matric, C , satisfies its own characteristic equation 

n 

l -0^=0 (a 0 =l) (5) 

1=0 

It follows from Eq, (5) that the entire n-space is annihilated. Thus, 

n 

^ ce C n L z ■ = 0 (6) 

i=b 1 


and 


(7) 


n 


V n~i 

) or. G z 

L— 1 



Z 


Equation (7) can be used to determine the n unknown coefficients, O'. . Once 
these are determined, we have, after multiplication by G ^ , 

c " lz = --til “i 0 " 1 " 2 } < 8 > 

n i=o 

To realize a saving in computer time, we require an approximate characteristic 
equation for the matrix C, of much lower degree than n„ Thus, we are led 
to find the p , 0 , , scalars which minimize the length of the vector, t. 

P 

£ J5. C 11 " 1 z = l , p« n (9} 

i=o 1 


The 0's are given by the least squares solution for -1=0, with /? 0 = 1, We 
form the Gram-Schmidt upper triangular decomposition of the nxp+1 matrix 

p 

of the vectors z, Cz, , . C z, We have 

(nxp+1) (nxp+1) 

(z , Cz, C 2 z , . . . , C P z ) = (0)(T) • (10) 

where (9) is an orthogonal matrix (nxp+1) and (T) is upper triangular 
matrix (p+lxp+1). 


The solution for 0 is given by 


(p+lxl) 


ft 


(t(p+1x P +1)) 


P 

J 

P-1 


= 0 


( 11 ) 


JS. 
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st* 

Transposing the p-KL vector of (T), we have, for ,6^ 


r J3 

p 

-1 

t l(p- ! 'l) 

0 i 

# 

* 

= “(t (pxp)) 

t 2(p-H) 

i 

(_! ■ . 


* 

- W+l) - 


( 12 ) • 


The first few solutions are provided below. 


Zero Order 


x = 


2 . 


good only if 




a.. 

31 


« 1 . 


First Order 


x 



z 


(13) 


£14) 


Second Order 


x 


t t 
23 11 


*23 ~ t 13 *22 


z + 


t ll t 22 


t t -t t 
12 23 13 22 


C z 


( 15 ) 


The Gram-Schinidt coefficients are listed below: 


t.. = 

T i-1 
0. C z 

13 

3 

t.. = 

C(c- z -y 

n 

3”1 
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4 11 

T 

(z 

Z) 

0. = 

_1_ 

i-1 

fc^.-T 

1 

t.. 

u 

l L 

3”1 


To obtain the upper triangular- inverse of T, let S be the inverse of T, 
then the nonzero upper triangular elements of S are given by 


i(t -.g 

k=l 


From the above, it follows that the £ coefficients are given by 


^p-KL-i S i,p+1 


i = 1,2, . „ . ,p 


Since we are interested only in the ratio “ , the t ^ coefficient need 

P 

not be generated and may be arbitrarily set equal to unity. 
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ABSTRACT 


A variable-metric algorithm is described that makes use of both linear and 
quadratic penalty terms for handling nonlinear constraints and employs both pro- 
jection and penalty features. Quadratic penalty coefficients are adjusted in a 
process which attempts to maintain a positive-definite matrix of second partial 
derivatives of the function including penalty terms without generating the large 
positive eigenvalues traditionally attending the use of quadratic penalties, which 
cause zigzagging and slowed convergence. The schemes contemplated make use 
of inferred second-order properties not only in terms of the variable metric of 
DFP (or its relatives) but by estimation of second directional derivatives by fitting 
cubics to various functions along directions of search. Some experiments are de- 
scribed with a simple constrained-minimum problem contrived to offer difficulties 
with methods that use only linear penalties, hence taxing the quadratic -penalty- 
adjustment procedure. 
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INTRODUCTION 


The arrival of variable-metric parameter optimization, the Davidon-Fleteher- 
Powell algorithm (Ref. 1) and its relatives, literally revolutionized numerical optimi- 
zation in the sixties. Even variational problems, crammed into the mold by some- 
times awkward parameterizations, were treated handily by DFP in competition with 
various sophisticated continuous -control algorithms. The key to success is the 
superficially first-order character of the technique — only first partial derivatives 
need be generated explicitly — together with speed and the sureness of convergence 
accomplished by inference of second-order properties. 

But in most variable-metric applications work the constraints are treated by 
the quadratic penalty function, a primitive device well known to affect convergence 
rate adversely and to magnify numerical errors. The combination of penalty function 
and variable metric was explored in a 1966 paper (Ref. 2) which Included various 
auxiliary devices to ameliorate the adverse effects of penalty-function approximation. 
This particular computational procedure has turned out to be a reliable work-horse 
and is currently in fairly wide use In day-to-day applications work. 

Efforts at adapting variable metrics to the standard alternative scheme for 
treating constraints, gradient projection, proved straightforward and immediately 
tractable only in the case of linear constraints (Ref. 3). Variable-metric projection 
schemes, making selective use of what amount to linear penalty functions, were 
eventually developed for the case of norm mar constraints and proved workable in 
limited tests (Refs. 4 and 5). This class of variable- metric scheme has only seen 
limited use in complex applications, however, and is not yet highly developed. 

One suspects that current-state-of-the-ait schemes are costly and slow 
compared to what is possible. The focus in the following is upon that class of 
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problems in which auxiliary vector-matrix computations are inexpensive in com- 
parison with the generation of function samples and gradients, as, for example, in 
aerospace trajectory-shaping problems. Use is made of both penalty and projection 
ideas and various other features of the algorithms of Refs. 1-10; the adjustment of 
penalty coefficients represents the main innovation, and the bulk of the discussion 
will be devoted to this, 

A VARIABLE-METRIC GRADIENT PROCESS WITH LINEAR-PLUS-QUADRATIC 
PENALTIES 


. Consider an alternative to the problem of minimizing a function f(x) 

{ x an n-vector) subject to an m-vector equality constraint g(x) - 0 , namely 
the minimization of the function f given by 

f = ftgU - g Sg T (1) 

which contains both linear and quadratic penalty terms. With X= 0 and the 
elements of the diagonal matrix k.. » 0 , one has the quadratic penalty scheme 
(Ref. 6); large k values are needed in this approach not only to insure that the 
function adopted for minimization actually possesses a minimum near the con- 
strained minimum sought, but also to render the magnitudes of the constraint 
.violations, |g. j , small at the minimum. 

Hestenes’ Method of Multipliers (Ref. 7) employs both linear and quadratic 
penalty terms, with the quadratic terms viewed as primary; the linear terms, 
missing in a first major iteration, are introduced as auxiliaries to reduce con- 
straint violations and permit use of somewhat lower quadratic penalty coefficients. 

rVJ 

The X vector for each major iteration, which consists of a minimization of f, is 
taken in this algorithm as the value of X+ Kg at the end of the preceding major 
iteration. Of course, any minimization algorithm can be used for the major itera- 
tions but, for such unconstrained problems, DFP^and its relatives are highly 
competitive. 
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■ The algorithm examined in the following makes use of the form of f above 

including both linear and quadratic penalty terms. However, the viewpoint taken 

is different from the Method of Multipliers, namely that the linear terms are 

primary, the quadratic ones supplementary and missing whenever advisable. The 

X- vector components will be determined as the projection values every few cycles, 

and the K diagonal elements chosen generally so as to provide the second partial 

derivative matrix f with positive definiteness, but without the excessive mai’gin 

xx 

traditionally furnished by large quadratic penalty terms, which hinders convergence. 
The linear penalty terms provide the means of reducing constraint violations to 
zero in case the constraints are compatible, 1. e. , the surfaces defined by g^ = 0 
have an intersection. 


The two algorithms examined employ DPP (Ref. 1) and its batch-processor 
DFP modification (Ref. 8} applied to f for major iterations. They bear a resem- 
blance to the Method of Multipliers, differing from It in the determination of X and 
k values. It is similar for the first few cycles, during which the diagonal K ele- 
ments are assigned "moderately large” positive values in the quadratic-penalty- 
function tradition. The major iteration proceeds by variable metric for n cycles, 
however, rather than all the way to a minimum. 


The general idea of the quadratic-penalty-coefficient selection scheme is 
control of the eigenvalues of the second-partial-derivative matrix 


f 

xx 


f + 

XX 


m 

I 

j=l 


g 3 




XX 


m 

I 




( 2 ) 


to produce positive-definiteness and a largest eigenvalue not much exceeding the largest 
eigenvalue of f + g ,X [illegal notation but suggestive shorthand for the first two 
terms of (2)1. One would like this not locally, with X the projection value, but at 
the constrained minimum where the projection X coincides with the Lagrange multi- 
'plier vector; however, it would be difficult and expensive enough to calculate the local 
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second partials and the largest eigenvalue, so a less direct and more approximate 
approach is taken. The scheme proposed as follows takes advantage of-the fact that 
there will usually be a large range of values for the k meeting the requirements, 
the lower limit determined by loss of definiteness and/or excessive constraint viola- 
tions, and the upper limit related to the largest eigenvalue of f + g X. 

XX XX 


During the n cycles of each "batch”, or major iteration, second directio nal 
derivatives along the n directions of search are estimated for the function 
f*.= f + g A* , where X* is given by 


-(SxVj 1 


T 

g H f 

X o X 


( 3 ) 


as the gradient -projection value; this varies from cycle to cycle. is a full-rank 
nxn matrix, fixed during a batch. At the constrained minimum sought, the value 
given by (3) is equal to the Lagrange multiplier vector for stationary f + g A; it is 
independent of the metric when evaluated at the constrained minimum. 

Tor a step determined by the modified DFP algorithm as 


Ax, = x - x. 

i l+l i 


i = 1, , n 


( 4 ) 


the first and second derivatives in the direction are given by 


f*' = 


. T 
Ax f 


x 


Ax 


(5) 


f* + 

f = "T&rr 


( 6 ) 


f* = 


6(f**- f*) _ 2f*' + +4f* t 
|Ax| 2 | Ax | 


(?) 


*”+ - 6 ( f* - f 8 ¥+ ) 2f*' + 4f i ' ,+ 


Ax 


Ax j 


(8) 
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(Here the + superscript denotes evaluation at the 1+1 end of a search segment. ) 
The second derivative estimate corresponds to cubic fit to f ’ and f ■ values at 
the endpoints. In computations carried out with short word length or subject to 
excessive round-off error, the simple difference-quotient approximation which is 
the average of (7) and (S) may be preferable. 

In the vicinity of the constrained minimum sought, some ox the second di- 
rectional derivatives of the function £* , which approximates the first two terms 

t** 

of f , can be expected to be positive as f +g X possesses at least n-m non- 

XX XX 

negative eigenvalues. The largest positive value determined over one or several 
batches can be adopted as a guide for adjusting the penalty coefficients, as it will 
fdll in the range between zero and the largest eigenvalue. 

The second directional derivatives of f* in directions along the constraint 
function gradients are not, in general, positive; if they were, in a large enough 
neighborhood of the constrained minimum, the quadratic penalty terms might be 
dispensed with. One set of requirements on the quadratic penalty coefficients 
might be determined from second derivatives of f in these directions, by requiring 
them to he equal at the least to a guideline value. 


Carrying this scheme out directly necessitates either special probing opera- 
tions in the directions of the constraint gradients or the inference of equivalent 
information from the function samples and gradients computed in the course of 
minimization iterations. Both have been considered and investigated in a preliminary 
way and a combination is recommended for use. An estimate of the latter type for the 
penalty coefficients k. is given by the maximum (over one or more batches) of the 
values k. given by 



i=l, , n 

3-1. » m 


( 9 ) 
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where. 



( 10 ) 


Values are to be excluded from consideration when the two terms of the denominator 
in (9) are opposite in sign and nearly equal in magnitude; likewise when (3 given by 
(10) is smaller than some prescribed value, indicating that the particular step Ax 
was nearly in the tangent plane of the constraint whose quadratic penalty coefficient 
requirement is being estimated. Here f = f 4- g X, where X is the value of the 
linear penalty X employed in the function f during the particular batch; a prime 
denotes the first derivative along the direction of the step taken, a double prime the 
second derivative. The expression (9) was obtained by requiring the k value be 
large enough to produce f” equal to the guideline value cf^_ for £=1 and re- 
main bounded for small j3 (inasmuch as the denominator behaves like for 
g - 0 and /3 small). Since it is desired that k estimates err on the high side, 
the f" values used should be the larger of the values at beginning and end of the 

• n “n 

search segment for f* and the smaller of the two values for f . 


An additional candidate is introduced to cover the frequently-occurring 


contingency that all are small over one or more batches used in the selection. 

3 i 

viz. 

, c *»i , 

(ei - f . ) 

. max mm 

k. = — 


3 


s j 


(ii) 


x 


* 

where & corresponds to the last cycle before k selection and f ' . is taken as the 

mm 

smallest of the f ' values over a chosen number of batches, or zero, whichever is 

the lesser. The multiplicative constant e a 1 in the guideline value of f ‘ introduces 

a measure of conservatism to offset the possibility that none of the candidate values of 

f * n in the maximization determining f is really close to the largest eigenvalue of 

f \* 

g xx^- * 
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METRIC ADJUSTMENT 


Aftex* determining the A. and K elements anew at the beginning of each 

major iteration, one would like to adjust the variable-metric matrix II to account, 

at least approximately, for the changes. The corrections are based upon the idea 

~ -1 

that the H matrix emerging from the preceding major iteration approximates f ^ . 

TOC 

No correction is made for changes in the sum A+ Kg appearing in the second 

partials (2) as this sum approximates the Lagrange multiplier at the constrained 

minimum when the g. are small. 

1 


Corrections for k. changes are done sequentially, using 


H + AH = 


W-) Hg i h 


1 + Alt. g. Hg. ' xx 

ii i 

x x 


which accounts for changes in the last term of eq. (2) via the Schur identity (Ref. 2). 
Each increment Ak. is limited to some fraction of the original or updated 

i 

value k. so as to insure that the denominator of the fraction in parenthesis re- 
mains positive and does not nearly vanish. 


TEST PROBLEM 


The problem used for experiments employed a cubic in one variable, x , 
for f, and a quartic of the following form for the single constraint function g : 


A 2 3 

X 1 +£L 1 X 1 +a 2 X l 


. 2 , 2,4 

X l' b l X 2 "Vs - b 3 X 2 
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•In the simplest ease, used for functional checks of computer programming, 
a - a = b = 0 , b 1 > 0 , b > 0 , the constraint surface is a paraboloid of elliptic 
cross-section and the minimum of the linear function f is attained at the origin. 

The constraint function nonlinearity is an essential feature of the well-defined con- 
strained minimum. If a slightly negative value of is introduced, one already 
has a problem for which no minimum of g A. exists at the constrained minimum 
as the Hessian matrix is indefinite and, accordingly, quadratic penalty terms are 
essential. This is not quite enough complexity for algorithm development, evaluation, 
and comparison, however, as f-rg X is then quadratic and the variable-metric pro- 
jection schemes have too easy a time of it. Hence, use of a f 0 and b ^ 0 is 

d tJ 

attractive. It should be noted that a g > 0 large enough precludes the appearance of 

minima other than at the origin. The numerical values of the coefficients used in the 

~2 “3 2 —1 

computational experiments were: a = -10 , a = 10 , b = 1 » b = 10 , b - 10 ; 

X d X St O 

these are such as to offer modest challenge. 

The starting point for the numerical computations of the example was x^ = 10 , 
x = 5 , and x = 10 . The multiplicative constant c used in (9) and (11) to designate 
the guideline value of r was taken as unity in the comparison. 


COMPUTATIONAL COMPARISON 

To afford a basis for comparison, DPP was run on the example with the quad- 
ratic penalty coefficient fixed at several values and zero linear penalty coefficient. 

The first three entries in the accompanying table present these results for quadratic 

3 2 

penalty coefficients of 10 , 10" 1 , and 10 , At the minima found, the constraint 
g=0 was not satisfied owing to the absence of linear penalty terms, 'boundary shifting' 
or any other palliative. The violations were found to be excessive for k ^ 10 , 
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The next two entries lli the table are for sequential and batch-processor 
versions of the algorithm described in the preceding. The batch version was un- 
accountably better than the sequential; usually, the sequential is slightly better 
with fixed quadratic penalty when an accurate linear search is performed (Ref. 8) . 
The accelerated search of Ref. 9 was employed in all of the computations presently 
reported, with a tight tolerance employed for termination. The quadratic-penalty 
adjustment procedure was restrained from changing the coefficient more than an 

order of magnitude on any single adjustment. The scheme brought the coefficient 

-X 

down into the range 10 ^ k ^ 10 , a favorable range when the linear term is 

present to avert large constraint violations. 

The last two entries correspond to the variable-metric projection algorithms 
of Refs. 5 and 4, respectively, the latter slightly modified. These are reviewed in 
Appendices A and B for the reader's convenience. 


CONVERGENCE COMPARISON 

Algorithm .... 

Quadratic Penalty 
Coefficient k 

Number of Cycles 
:-to Convergence 

DFP 

io 3 

105 

DFP 

io 2 

58 

DFP 

10 

24 

linear-quadratic penalty/sequential 

variable 

27 

linear-quadratic penalty /batch 

variable 

21 

Rosen-ICreuser (modified) 

projection 

64 

Kelley-Speyer 

projection 

72 
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POSSIBLE IMPROVEMENTS IN PENALTY-COEFFICIENT-ADJUSTMENT 


PROCEDURES 


Two features intended to aid the process of adjusting linear and quadratic 
penalty coefficients •.=. ill be described very briefly. These have been explored 
computationally and found to produce reasonable results, but evaluated inyufuoiently 
to permit overall judgment on their merits. 

Early values of the linear penalty coefficients, the components of by 
projection, tend to be far off the mark, which recommends use of a zero value 
during the first batch. A short first batch suggests itself, m cycles instead of n , 
mainly to reduce the magnitudes of the constraint violations. It might be hoped that 
a better metric would also be obtained as well, better at least in the subspace de- 
fined uy the constraint gradient vectors. An obvious modification of the batch metric 
update formula of Ref. S is required. 

Another obvious temptation is smoothing of X values used on successive 
batches by weighted averaging, heavily weighting the new projection value when 
there is an indication of accelerating convergence, as by drastic shrinkage in the 
magnitude of the projected gradient vector during the batch just completed. The 
motivation for avoiding unduly large fluctuations in linear penalty coefficients, of 
course, is that changing the function being minimized taxes the machinery for 
inferring the metric and the various second derivatives. 

In the limited trials of these two features to date, it has been found that 
they generally enhance the smoothness and "surefootedness" of the algorithm, al- 
though at slight expense in convergence speed. The use of higher values of the 
constant c a 1 , say 2 , or even 10 , has a similar effect. 
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CONCLUDING REMARKS 


The results are thought to indicate promise for the class of algorithm com- 
bining linear and quadratic penalty adjustment with variable -metric optimization. 
More extensive testing is obviously needed, including the large class of problems 
in which the quadratic terms can safely be adjusted downward to zero. 
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APPENDIX A 


THE ROSEN-KREUSER PROJECTION ALGORITHM 


After restoration of constraints, the gradients of f and g and the projection 

a „ A 

multiplier vector are calculated and designated f , g , X at this batch-reference 

X X 

point x = x. 


a m —I A rn a 

X = - (g H g ) g H f 

K X O X X O X 


(Al) 


The algorithm proceeds to minimize f + g X subject to the linear constraint 


gjx-x) = 0 

A 


{A 2) 


by projecting the gradient of f + g X upon this constraint at each step. The pro- 
jection multiplier needed is 

aT a. ”1- aT * 

X = - (g H g ) g X H (f +g X) (A 3) 

X o X O V X X 

A linear search is made in the direction 


Ax = - O' H ( f x 4 'S x X+g x X) 


(A 4) 


to a one-dimensional minimum. On the first cycle X= 0 , but not subsequently, 
except in special cases such as linearly constrained problems. The metric H is 

A 

updated sequentially by the DEP formula evaluated using the gradient of f + g X: 


H(Af + Ag X)(Af : + Ag X) T H 
H + AH = H - — — — + 


Ax Ax 


(Af x + Ag ;( X) T H (Af x + Ag x X) Ax T (A f x + Ag x X) 


(A 5) 


After n-m cycles, H attains its limiting value for a quadratic f, linear 
g model; hence n-m is a natural batch size. It would seem generally more ef- 
ficient to restore and relinearize after each n-m cycles of DFP than to run to a 

A 

minimum of f + gX as proposed in Ref. 5. In fact, this feature encountered diffi- 
culty in the numerical computations reported, and relinearization each n-m 
cycles was the modification actually used. 
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APPENDIX B 


THE KELLEY-SPEYER PROJECTION ALGORITHM 


The accelerated gradient projection process of Ref. 4 employs the formulas 


Ax = - a H ( f x +£ x *) (B 1) 

X - -(B^Hgp ^Ef x (B2) 

\vith X recalculated every optimization cycle which successfully terminates on 
a one-dimensional minimum of f + g X , and with H updated by 


H + AH = H -i- 


T 

Ax Ax 

Ax T (Af x +Ag x A) 


H(Af x + Ag^ X)( Af + Ag^ X) T H 
(Af +Ag X) T H(Af +Ag X) 

A X A A 


(B 3) 


which is the DFP formula applied to the linear combination f + gX, hence 
guarantees that I-I remains positive definite. Constraint restorations are 
carried out after each optimization cycle (Ref. 10). 
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